
    GetName = function(var) deparse(substitute(var))
    GetVar  = function(varname)  eval(as.symbol(varname))
    
    logit = function(p) log(p/(1-p))
    expit = function(x){
        output = rep(1,length(x))
        index  = !is.na(x) & x <= 100
        e.x    = exp(x[index])
        output[index] = e.x/(1+e.x)
        output[is.na(x)] = NA
        return(output)
    } 
    
    Mean = function(x) mean(x, na.rm=TRUE)

    NameVar = function(list,name){
        
        varname1 = strsplit(name,split=".",fixed=TRUE)[[1]][1]
        paste(varname1,"=",list[name][[1]],sep="")
        
    }
    
    NameList = function(x){
        
        sapply(1:length(x), function(i) NameVar(x,names(x)[i]))
        
    }
    
    NameFirst = function(names){
        
        sapply(names,  function(name) strsplit(name,split=".",fixed=TRUE)[[1]][1])
        
    }
    
    ParaValues = function(paras){
        
        ###     Create an array containing all the parameter value combinations
        para.string = paras[[1]]
        for(i in 2:length(paras))   para.string = outer(para.string,paras[[i]],"paste")
        
        ###     Extract numerical values from strings
        Extract = function(name) as.numeric(strsplit(name,split=" ",fixed=TRUE)[[1]])
        para.value = apply(para.string,1:length(dim(para.string)),Extract)          # The first dimension contains the results
        dimnames(para.value)[[1]] = NameFirst(names(paras))
        dim.order  = c(2:length(dim(para.value)),1)
        return(aperm(para.value, dim.order))
        
    }
    
    MyAttach = function(list){  # attach a list to the global environment
        for(i in 1:length(names(list))){
            assign(names(list)[i], list[[i]], envir=.GlobalEnv)
        }
    }
    
    
    Message = function(paras, para, ptm1){
        para.min = sum(sapply(paras[1:(length(para)-1)], min))
        output   = NULL
        if(sum(para[1:(length(para)-1)]) == para.min){
            if(para["n"] > min(paras$n.cand)){
                output = paste(output, " takes ",round((proc.time() - ptm1)[3],
                                                       3),"s \n",sep = "")  
                ptm1 <= proc.time()
            }
            output = paste(output, "Sample size ", para["n"], sep = "")
        }
        return(output)
    }
    
    # Print = function()
    
   
    
    
    ##### Rounding functions (start) #####
    
    RoundS = function(x,num){
        #' Round a number by scientific digits
        #'
        #' Args:
        #'     x: The number to be rounded
        #'     num: the number of scientific digits
        #' Returns:
        #'      The rounded number
        #'
        gsub("\\.$", "", formatC(signif(x,digits=num),digits=num,
                                 format="fg", flag="#"))
    }
    
    RoundD = function(x,num){
        sprintf(paste("%.",num,"f",sep=""), round(x,num))
    }
    
    Round = function(x, num){
        # Round large numbers by scientific digits
        # Round small number by decimal places
        output = Round(x,num)
        output[x<1 & !is.na(x)] = Round_decimal(x[x<1 & !is.na(x)],num)
        return(output)
    }
    
    ##### Rounding functions (end) #####
    
    
    
    Optimize = function(startpars,objective,Hessian = FALSE){
        # input
        # startpars:          parameters to start with
        # neg.log.likelihood: objective function to minimize
        # Hessian:            logical; whether or not calculate the Hessian
        
        opt = optim(startpars,objective,hessian=Hessian)
        if(opt$convergence != 0){
            if(MESSAGE) cat("The default method does not converge in 500 steps. \n")
            opt = optim(startpars,objective,control=list(maxit=10000))
            if(opt$convergence != 0){
                if(MESSAGE) cat("The default method does not converge in 10,000 steps... \n")
                opt = optim(startpars,objective,control=list(maxit=100000))
                if(opt$convergence != 0){
                    if(MESSAGE) cat("The default method does not converge in 100,000 steps.. \n")
                    if(MESSAGE) PrintErrorMessage(opt)
                }else{
                    if(MESSAGE) cat("But it converges in",opt$counts[1], "steps...\n")
                }
            }else{
                if(MESSAGE) cat("But it converges in", opt$counts[1], "steps...\n")
            }
        }    
        
        return(opt)
    }
    
    
    ### 
    
    GetfZ = function(f1, f0, Z){
        fZ       = f1
        fZ[Z==0] = f0[Z==0]
        return(fZ)
    }
    
    
    getProb = function(theta, phi1, phi2, phi3, phi4, op, target){
      
      if(target == "LATE")  f0f1 = mapply(getProbScalarRD,atanh(theta), log(op))
      if(target == "MLATE") f0f1 = mapply(getProbScalarRR,log(theta), log(op))
      if(length(op)>1) {   f0   = f0f1[1,]; f1 = f0f1[2,]  }
      if(length(op)==1){   f0   = f0f1[1];  f1 = f0f1[2]   }
      
      p011 = (1-phi1) * (1-phi2) * phi3
      p001 = (1-phi1) * (1-phi2) * (1-phi3)
      p110 = (1-phi1) * phi2 * phi4
      p100 = (1-phi1) * phi2 * (1-phi4)
      p111 = f1 * phi1 + p110
      p010 = f0 * phi1 + p011
      p101 = 1 - p001 - p011 - p111
      p000 = 1 - p010 - p100 - p110
      
      if(length(op)>1)  return(cbind(p000,p001,p010,p011,p100,p101,p110,p111))
      if(length(op)==1) return(c(p000,p001,p010,p011,p100,p101,p110,p111))
      # if(length(op)==1) return(list(p000=p000,p001=p001,p010=p010,p011=p011,
      #                               p100=p100,p101=p101,p110=p110,p111=p111))
    }